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Abstract 


In this paper we develop two goal-oriented adaptive strategies for a posteriori error estimation within 
the generalized multiscale finite element framework. In this methodology, one seeks to determine the 
number of multiscale basis functions adaptively for each coarse region to efficiently reduce the error in 
the goal functional. Our first error estimator uses a residual based strategy where local indicators on each 
coarse neighborhood are the product of local indicators for the primal and dual problems, respectively. 
In the second approach, viewed as the multiscale extension of the dual weighted residual method (DWR), 
the error indicators are computed as the pairing of the local H~^ residual of the primal problem weighed 
by a projection into the primal space of the Hq dual solution from an enriched space, over each coarse 


neighborhood. In both of these strategies, the goal-oriented indicators are then used in place of a standard 
residual-based indicator to mark coarse neighborhoods of the mesh for further enrichment in the form 
of additional multiscale basis functions. The method is demonstrated on high-contrast problems with 
heterogeneous multiscale coefficients, and is seen to outperform the standard residual based strategy with 
respect to efficient reduction of error in the goal function. 

1 Introduction 

Many practical problems are multiscale in nature, including flow in porous media, seismic wave propagation, 
and physical processes in perforated media. These problems are described by partial differential equations 
(PDEs) with potentially high contrast multiscale coefficients. Direct computation of high resolution discrete 
solutions to these problems can be very expensive. Typically, some type of model reduction techniques 
are used to solve multiscale problems. Established techniques include numerical homogenization methods 
[mini 12] and multiscale methods miiMiiiiiiiiiiMiisiiiiiiiTiiiiiiiin]. In numerical homogenization 
methods, the upscaled media properties are computed over coarse grid blocks, each of which is much larger 
than a characteristic length scale. In inultiscale methods, local multiscale basis functions determined by 
local fine-scale problems are constructed in each element of the coarse grid. Generally, one uses a few 
multiscale basis functions in each coarse element to approximate the global solution by solving a coarse 
mesh problem over the entire domain. Multiscale basis functions are constructed in an offline step before 
the coarse mesh problem is solved, after which some type of adaptivity is needed to choose multiscale basis 
function appropriately. 

In this paper we will use multiscale methods, where multiscale basis functions are constructed in each 
coarse region, as illustrated schematically in FigurejU To be more specific, we consider a multiscale problem 


L{u) = /, 

where L = —div(A^(x)V'u) and seek the solution in the form 



* Department of Mathematics, The Chinese University of Hong Kong, Hong Kong SAR. Eric Chung’s research is supported 
by the Hong Kong RGC General Research Fund project 400813. 

^Department of Mathematics, Texas A&M University, College Station, TX 
^Department of Mathematics, Texas A&:M University, College Station, TX 


1 



In each coarse region Wj, we construct a set of multiscale basis functions , i = These multiscale 

basis functions, as described below, will be constructed in the offline stage using the generalized multiscale 
finite element method (GMsFEM) and represent the local heterogeneities of the solution space. In earlier 
works on the multiscale finite element method (MsFEM) [35], the authors sought one multiscale basis function 
per coarse element. However as it was later argued, one may need additional basis functions over each coarse 
element for a sufficiently high hdelity approximation. It is further shown that the optimal number of 
multiscale basis functions in each region depends on the heterogeneities in the solution space. Typically, 
some adaptive criteria based on a posteriori error estimation is used to determine how many basis functions 
to choose in each coarse region uii. For exampe, in m. the authors develop an error indicator based 
on the norm of the residual to determine the number of basis functions to add in each region over 
each adaptive iteration. Adaptive multiscale methods follow traditional adaptivity concepts as in miEiiTi 
|3Zl ISOjlMl io]; however, the multiscale indicators contain information about local heterogeneities. Earlier 
approaches to goal-oriented adaptive methods for multiscale problems include numerical regularization and 
numerical homogenization with adaptive mesh refinement as in [381134] . To the authors’ knowledge, the 
current presentation is the first to develop a goal-oriented enrichment strategy within the general GMsFEM 
framework. 

For many practical problems, one is interested in approximating some function of the solution, known 
as the quantity of interest, rather than the solution itself. Examples include an average or weighted average 
of the solution over a particular subdomain, or some localized solution response. In these cases, goal- 
oriented adaptive methods yield a more efficient approximation than standard adaptivity, as the enrichment 
of degrees of freedom is focused on the local improvement of the quantity of interest rather than across the 
entire solution [Ml II IMl m mi [Ml H ISl E] • In this paper, we study goal-oriented adaptivity for multiscale 
methods, and in particular the design of error indicators to drive the adaptive enrichment based on the 
goal function. In multiscale methods, goal-oriented adaptivity can play an important role in the efficient 
approximation of the quantity of interest as heterogeneities in the coefficients may require standard adaptive 
methods to add degrees of freedom in regions with limited influence on the goal function. In this paper, 
we develop a goal-oriented approach for multiscale methods within GMsFEM framework. In the proposed 
approach, we increase the accuracy of the approximation by enriching the space rather than refining the 
mesh by choosing multiscale basis functions computed in the offline stage. 

For multiscale basis construction, we use GMsFEM. The construction of multiscale basis functions uses 
local snapshot spaces and requires solving local spectral problems over each coarse element. The local 
snapshot functions represents the solution space in each coarse region, and they can include all possible local 
fine-grid functions or harmonic functions. In the snapshot space, we perform a local spectral decomposition 
and select multiscale basis functions which correspond to the dominant eigenvalues. The multiscale basis 
functions are constructed by multiplying the dominant eigenmodes by a partition of unity function, e.g., 
multiscale partition of unity function. In [16] , we developed an adaptive approach and developed a posteriori 
error indicators, which include the information from local spectral problems, e.g., the value of the eigenvalue 
corresponding to the first eigenvector not included in the coarse space. We derived error estimates and 
presented numerical results which demonstrate the improved efficiency of the adaptive approach, guided 
these indicators. For goal-oriented problems, we now design goal-oriented error indicators, which are different 
from those developed earlier m for multiscale problems, by the additional consideration of a dual problem 
to direct the adaptivity towards the approximation of the quantity of interest. 

In this paper we develop two goal-oriented adaptive strategies for a posteriori error estimation. Our first 
error estimator uses an idea similar to a standard residual based adaptive method, and can be seen as the 
multiscale extension of the fip-adaptive method presented in [^. In this case the elementwise indicator is 
formed by the product of local residual indicators for the primal and dual problems, respectively. In the 
second approach, viewed as the multiscale extenstion of the dual weighted residual method (DWR), the 
error indicators are computed as the pairing of the local H~^ residual of the primal problem weighed by a 
projection into the primal space of the dual solution from an enriched space. In both of these strategies, 
the goal-oriented indicators are then used in place of a standard residual-based indicator to mark coarse 
elements of the mesh for further enrichment in the form of additional multiscale basis functions. 
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The remainder of the paper is organized as follows. In Section [5] we present and overview of multiscale 
methods and adaptivity. In Section [3] we give a detailed description of the construction of the multiscale 
basis functions, review some results on residual-based adaptivity for GMSFEM. In Section U we introduce 
two goal-oriented a posteriori error indicators and present an algorithm for goal-oriented adaptivity. Finally 
in Section [S] we present numerical results demonstrating the efficiency of the proposed goal-oriented error 
indicators. 


2 Overview of Concepts 

In this paper, we consider second order multiscale elliptic problems of the form 

—div(K(a;) Vu) = / in D, 
u = 0 on dD, 


where D is the computational domain, and k(x) is a scalar valued heterogeneous coefficient with multiple 
scales and high contrast. The problem m can be solved by many classical numerical techniques, such as the 
conforming finite element method, but with extremely high computational complexity due to the fact that 
a very fine mesh is necessary to resolve the multiscale nature of the solution. Thus, some multiscale model 
reductions are needed to compute an accurate solution efficiently. In the following, we give a brief overview 
of GMsFEM and its basis enrichment techniques as applied to problem ©■ Let u GV = Hq{D,) be the true 
solution satisfying 

aiu,v) = if,v), vGV, (2) 

where a{u, v)= k{x)'Vu ■ Vv dx, and (/, v) = I fv dx. Define the energy norm on V by ||u||y = a{u, u). 

J D J D 

To introduce the GMsFEM for the problem m, we first give the notion of fine and coarse grids. We let 
be a standard conforming triangulation of the computational domain D into finite elements, which can 
be triangular, rectangular or some other polygons. We refer to this partition as the coarse grid. Subordinate 
to the coarse grid, we define the fine grid partition, denoted by by refining each coarse element into 
a connected union of fine grid blocks. We assume the above refinement is performed such that 7”^ is a 
conforming partition of D. We let N be the number of interior coarse grid nodes, and let {sijfci be the set 
of coarse grid nodes or vertices of the coarse mesh ■ Moreover, we define the coarse neighborhood of the 
node Xi by 

uJz = e T^; Xi e Kj}, (3) 

which is the union of all coarse elements which have the node Xi as a vertex. See Figure [T] for an illustration 
of the coarse elements and coarse neighborhoods within the coarse grid. We emphasize the use of u)i to 
denote a coarse neighborhood, and K to denote a coarse element throughout the paper. 

Next, we briefly overview the continuous Galerkin (GG) formulation of GMsFEM, a generalization of 
the classical MsFEM [32]. For each coarse node we define a set of basis functions supported on the 
neighborhood Wi. We denote the fc-th basis function supported on the coarse neighborhood uJt by '0^% We 
remark that in the GMsFEM, we will use multiple basis functions per coarse neighborhood, and the index 
k represents the local numbering of these basis functions. These multiscale basis functions are constructed 
from a local snapshot space and a local spectral decomposition defined on that snapshot space. The snapshot 
space contains a collection of many basis functions that can be used to capture most of the fine features 
of the solution, and the multiscale basis functions are constructed by selecting the dominant modes 
of a local spectral problem. Using the these multiscale basis functions, the GG solution is represented as 
Ums{x) = Once the basis functions are identified, the GG global coupling is given through 

the variational form 

a(ui„s, v) = (/, v), for all v G Kfr, (4) 

where VqB is the space spanned by the basis functions and a(-, •) is the usual bilinear form correspond¬ 

ing to O- We remark that one can use other formulations, such as the discontinuous Galerkin formulation 
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Figure 1: Illustration of a coarse neighborhood and a coarse element 

(see e.g., [121 El ED]), the mixed formulation (see e.g., [131 E]) or the hybridized discontinuous Galerkin 
formulation (see e.g., [35]) to couple the multiscale basis functions. 

In using GMsFEM, it is desirable to determine the number of basis functions per coarse neighhorbood 
adaptively based on the heterogeneities of the coefficient k,(x) in order to obtain an efficient representation of 
the solution. In [16] , a residual based a posteriori error indicator is derived and an adaptive basis enrichment 
algorithm is developed under the CG formulation. In particular, it is shown that 

N 

IIU Mmsil V — ^ ^ ) G I 

i=l 

where is the residual of the solution Ums on the coarse neighborhood uii. Thus, local residuals of the 
multiscale solution give indicators to the error of the solution in the energy norm, and one can add basis 
functions to the coarse neighborhoods when the residuals Vi are large. Convergence of this adaptive basis 
enrichment algorithm is also shown in m- On the other hand, for some applications one needs to adaptively 
construct new basis functions in the online stage in order to capture distant effects. In m, such online 
adaptivity is proposed and mathematically analyzed. More precisely, when the local residual ri is large, one 
can construct a new basis function (j) € Vo{uJi) in the online stage by solving 

a{(l),v) = {ri,v), VvGVoiuJi), 

where Vo(wi) is the restriction of V in uJi with zero trace on dcui. Numerical results in [15] show that a couple 
of these online basis functions can help to reduce the error by a large amount. 

The adaptivity procedures discussed above are designed with the aim of reducing the error in the energy 
norm. In some applications, one may be more interested in reducing error measured by some function of the 
solution other than a norm. For example, in flow applications, one needs to obtain a good approximation of 
the pressure in locations where the wells are situated. Therefore, we now consider goal-oriented adaptivity 
within GMsFEM. Specifically, we define a linear functional g : F —>■ R, referred to as the goal functional. 
In goal-oriented adaptivity, one wants to adaptively enrich the approximation space in order to reduce the 
goal error defined by g{u — rtms)- In the construction of goal-oriented adaptivity for GMsFEM, we use local 
indicators based on the solution of a dual problem: finding z € V such that 

a* {z,v) = g{v), \/vGV. (5) 

For a primal problem a{u,v) = (/, u) based on bilinear form a( •, •): 1^6 dual form a*( •, •) is the formal 
adjoint of the primal, satisfying a*(w,v) = a(v,w), and in the current symmetric case, a*( •, •) is identical 
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to the primal. Formally, the primal-dual equivalence follows for u the solution to the primal problem (U) 
and z the solution to the dual problem ([5]) 

/(z) = a(u,z) = a*{z,u) = a(z,u) = g{u). (6) 

Error estimates for the quantity of interest g{u) follow from (jH]) and Galerkin orthogonality with respect to 
the discrete problems and their respective solutions. Forming error indicators based on both primal and dual 
problems and these estimates, we add multiscale basis functions to coarse neighborhoods when the values 
of the corresponding indicators are large. Our numerical examples show that the goal-oriented approach 
performs better than the residual approach for high-contrast problems when the error is measured by the 
goal-functional, g(u — Ums)- 


3 The GMsFEM and residual-based adaptivity 

In this section, we will give a detailed description of the GMsFEM (see for example [T9l[21]) and it’s residual 
based adaptivity (see for example M)- 


3.1 Local basis functions 


We first present the construction of the multiscale basis functions. This construction is performed in the 
offline stage; that is, basis functions are pre-computed before the actual solve of the problem. The construc¬ 
tion starts with a snapshot space. This space contains a relatively large set of basis functions which can 
be used to capture most features of the hne-scale solution. The next step is to perform a local dimension 
reduction to obtain a lower dimensional subspace that can still be used to approximate the solution with 
good accuracy. The local dimension reduction is performed by solving a spectral problem, and the dominant 
eigenfunctions are used as the multiscale basis functions. 

First, we define a snapshot space where the functions in L;n*ap are supported in Wi. The snapshot 

space can be the space of all fine-scale basis functions V{uji) = {v\^^ \ v &V} or the solutions of some local 
problems with various choices of boundary conditions. For example, we can use the following K-harmonic 
extensions to form a snapshot space. Specifically, let {i}}, j = --Li, index the set of fine-grid vertices 

that he on the boundary of each coarse neighborhood, duii. Define the unit source functions 5^(x) = 5(x®) 
for each j = 1,... Then construct the snapshot function g V{u]i) by solving 

—div(K(a;)V'0“‘’'*”‘‘'’) = 0, in wg 

= S!}, on du;,. 

The snapshot space corresponding to the region wg then contains Li functions 

y--,p = span{V>7-^"^P: 1 < J < L,}- 

We define the corresponding change of variable matrix 


TDl 

"^snap 


[ip 


uji ,snap 
1 5 




where are considered as the columns of the matrix. 

We next determine a set of dominant modes from L,nlip, and the resulting lower dimensional space is 
called the offline space To construct the offline space , we perform a dimension reduction of the 

space of snapshots using an auxiliary spectral decomposition. The analysis in |23| motivates the following 
generalized eigenvalue problem for eigenvalues Xjf and eigenfunctions in the space of snapshots: 




( 8 ) 


where 


/ 

J UJ 


= KZ] = / = (Kn^yARl 


snap 5 
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and 


= [S±] = [ = (i?Lap)^^i?Lap, 

J UJ 

where A and S denote analogous fine-scale stiffness and mass matrices as defined by 

Aij = / K.{x)V4>i ■ V(f>j dx, Sij = / 'R{x)<j)i(f)j dx, 

J D J D 

where (j)i is the fine-scale basis function for V. We will give the definition of k(x) later on. To generate 
the offline space we then select the smallest li eigenvalues from Equation ([5]) and form the corresponding 
eigenvectors in the space of snapshots by setting (for k = 1,..., where \E'^® 

are the coordinates of the vector and li is the number of eigenvectors chosen to span the offline space. 
We will use the set of local basis functions to form the approximation space in the next section. 

3.2 CG formulation 

In this section we create an appropriate solution space and the variational formulation for a continuous 
Galerkin approximation of Equation ©■ The idea is to use the basis set * = 1, 2, • • • , to form 

the approximation space, called the offline space, and apply the standard continuous Galerkin formulation. 
We begin with an initial coarse space 10 “'* = span{xi}^j^, where we recall N denotes the number of coarse 
neighborhoods corresponding to interior coarse nodes. Here, Xi 8-re the standard multiscale partition of unity 
functions which are supported in uji and are defined by 

-div (^(a;) Vxi) = 0, in K C oji, (9) 

Xi= Qi, on dK\duJi, 

Xi = 0, on duji, 

for all coarse elements K C uJi, where gi is a continuous function on dK which is linear on each edge of dK. 
Based on the analysis in |23j . the summed, pointwise energy k required for the eigenvalue problems (|8]) is 
defined as 

N 

i=l 

where H denotes the coarse mesh size. 

The partition of unity functions Xi 8re then multiplied by the eigenfunctions {4’k'’°^}k=i construct 
the multiscale basis functions 


i’i,k = Xi'^k'’°^^ 1 < * < and 1 < k <li, (10) 

where we recall k denotes the number of offline eigenvectors that are chosen for each coarse node i. We 
note the construction in Equation OT yields continuous basis functions due to the multiplication of offline 
eigenvectors with the initial (continuous) partition of unity. Next, we define the continuous Galerkin spectral 
multiscale space as 

VoS = spanji/ji^fc : I < i < N and 1 < fc < k}. (11) 

Using this offline space, we obtain the GMsFEM as in (jT]). 

3.3 Residual-based adaptivity 

In this section, we give a brief review of the adaptive basis enrichment algorithm proposed in |16j . After 
solving the coarse mesh problem and computing error indicators, the standard Dorfler marking strategy is 
applied with respect to the neighborhoods indexed by vertices i = 1,..., A^, rather than the elements K. 
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The marked neighborhoods are then enriched with additional basis functions. The spectral problem (|S|) gives 
a natural ordering of basis functions in the local snapshot space with respect to the eigenvalues, in 

increasing order. The analysis of m then suggests an adaptive procedure to add basis functions based on 
a local error indicator. The error indicator is defined by a i?~^-norm based residual, which gives a robust 
indicator with good performance for cases with high contrast media. 

Let uji be a coarse neighborhood and write Vi = Vo{u!i). For a given multiscale solution Ums, we define a 
linear functional i?“(u) on Vi by 


RUv) 



AcV-Ums • Vu, 


v€V. 


( 12 ) 


The norm of is defined as 


WR'^Wv- = sup 

vGVi 


Iklk. ■ 


(13) 


where ll^lly. = L K(x)|Vupc?x. The norm II 11 V* gives a measure on how well the solution u^s satisfies 
the variational problem (I2|) restricted to Vi. The norm ||i?“||y. can be obtained by solving an auxiliary 
problem of finding w G Vi such that a{w,v) = Rf{v) for all u S 1^ and then setting ||i?“||y = ||w||y. One 
can solve this auxiliary problem in the snapshot space to reduce the computational cost. Moreover, in 
m the following a posteriori error bound is proved 


N 

Ik-«n.sfy < CerrE + (14) 

where Cerr is a uniform constant independent of the contrast of n, and is the j-th eigenvalue of the 
eigenvalue problem (|8]) for the coarse neighborhood uJi. In particular, is the first, i.e., smallest eigenvalue 
from the spectral problem (jH]) for which the corresponding eigenvector is not included in the construction of 
the offline space. 

Using the above error bound (I14L a convergent adaptive enrichment algorithm is developed in [16) . We 
now describe this algorithm. The algorithm is an iterative process, and basis functions are added in each 
iteration/level based on the magnitudes of local residuals. We use m > 1 to index the enrichment level 
and let V^ be the solution space at the m-th iteration. For each coarse region, let be the number of 
eigenfunctions used at the enrichment level m for the coarse region cui. 

3.4 Adaptive algorithm 

The adaptive enrichment algorithm for GMsFEM is now summarized below. Choose a fixed marking pa¬ 
rameter 0 < 0 < 1. Choose also an initial offline space V^g by specifying a fixed number of basis functions 
for each coarse neighborhood, and this number is denoted by l], for each i = 1,... ,iV. Then, generate a 
sequence of spaces V^ and a sequence of multiscale solutions obtained by solving (g]). Specifically, for 
each m = 1, 2, • • •, perform the following calculations: 

Step 1: Find the multiscale solution in the current space. That is, find G V^ such that 

o-{u^s,v) = {f,v) foralluGUJg. (15) 

Step 2: Compute the local residual. For each coarse region a;,, compute 

,7? = iii?rii^,.(ArA+i)-\ 

where 

RUv) = [ fv- [ n{x)Vu^, ■ Vv, 

J UJi J U)i 


1 



consistent with m, and the norm is defined in m respectively. Next, re-enumerate the coarse 
neighborhoods so the above local residuals rjf are arranged in decreasing order r]f > r?! > • • • > Vn- 
That is, in the new enumeration, the coarse neighborhood wi has the largest residual rj^ and the coarse 
neighborhood ujn has the least residual 77 ^. 

Remark 3.1. An alternate approach to avoid the log TV complexity of the full sort is the standard 
binning or heapifying strategy Let rf = consider only rg that satisfy rj^ > (1 — 

6)%^/N. Let M = maxj rj^, and perform a partial sort of the remaining indicators collecting or binning 
the indices for which 2 ~pM < nf < for p = 0, 2 ,... ( 7 , where q is the smallest integer to 

satisfy 2 -('?+i)M < (1 - e)g^/N. 

Step 3: Find the coarse regions where enrichment is needed. Choose the smallest integer k such that 

N k 

i=l i=l 

The coarse neighborhoods a;i,a; 2 ,-'- , 0 ;*,, are then enriched with additional basis functions. If the 
partial sort of Remark 13.II is used in place of the sort, elements are marked by emptying the first bin, 
those indicators with M < rff < M/2, and then continuing on to the second bin, and so forth until (ITSl) 
is satisfied. As elements within bins are not sorted, this yields a quasi-optimal marked set, i.e., the 
marked set may not be the set of least-cardinality to satisfy (USD as in the full sort, but it is within a 
factor of two of the least cardinality. 

Step 4: Enrich the space. For each i = l,2,---,fc, add basis function for the region iVi according to the 
following rule. Let s be the positive integer such that Ajm+s-i-i is large enough compared with 
(see Remark . Then include the eigenfunctions • • • , 'L°m+s in the construction of the basis 

functions. The resulting space is denoted as . Mathematically, the space is defined as 

+ spanuti 

where nnd , with j = T™ -|- 1, • • ■ ,1/^ + s, denote the new basis 

functions obtained by the eigenfunctions • • • , In addition, we set + *■ 

Remark 3.2. The mathematical analysis in m specifies the choice of s. In practice, one can take 
s = I since the eigenvalues in m have fast growth. 

4 Goal-oriented adaptivity 

In this section, we present a goal-oriented adaptive enrichment algorithm for GMsFEM. The goal-oriented 
variant of the adaptive method requires the solution of a dual problem in addition to the primal at each 
iteration. The indicators are computed with both the primal residual and either the dual residual or a 
projection of an enriched dual solution into the primal space. These indicators predict which neighborhoods 
to enrich to increase the quality of the approximation of the quantity of interest. After introduction of the 
discrete dual problem, the finite dimensional analogue of ([5]), we propose two error indicators for goal-oriented 
enrichment. 

The dual problem plays a vital role in goal-oriented adaptivity as the vehicle for introducing the goal 
functional g into the adaptive process. Given a goal functional (/ : V —)> R, we define the discrete dual 
problem on approximation space VoS C V as: find z € I4ff such that 

a{v, z) = g{v), Vu S Foff- (17) 

As in ®, the discrete dual form a*( •, •) is identical to the primal a( •, •) for symmetric problems. The 
dual problem, however, features the goal functional g as the source. The discrete dual solution, may now be 
used to define goal-oriented error indicators. 


H ^-based goal-oriented indicator 

Our first goal-oriented indicator is similar in form to the residual based indicator described in the previous 
section. To motivate this indicator, we introduce the local bilinear form a{u,v)i = /t(x)Vu • ^vdx, and 
the induced localized energy norm ||u||yj = a{v,v)i. Let u be the solution to ([ 2 |), Ums G Kjff be the solution 
to (m, z the solution to JS]), and Zms G VoS the solution to (fT71) . Using Galerkin orthogonality and the 
relation between the primal and dual problems, the error in the quantity of interest satisfies 

g{u - Ums) = a{z, u - Ums) = a{u - Ums, z) = a{u - Ums, Z - Zms)- ( 18 ) 

Decomposing the global integration into neighborhoods by the partition of unity functions Xi given by ([9|) 

n „ n 

a{u-Ums,Z- Zms) ='Y' / - Ums) ' ^ {z - Zms) dx < \\u - UmsWy ^ \\z - ZmsWy ^ ( 19 ) 

i=l i=l 

where we used the fact that Xi is supported in uJi and the fact that |xi| < 1. Instead of using the norm of 
local residual for the primal problem defined in (fOl) to be our indicator, (fT91) suggests using the product of 
norms of local residuals for the primal and dual problems, posed in the same discrete space, Uoff- The local 
dual residual i?f : U —>■ R is defined by 


RKu) = 9{v) - K{x)Vzms 

OJi. 


where Zms S VoS is the solution to (IT71) . Analogous to (fOl) , the H ^ norm of Rf is defined as 

\\di^ \\V‘ - sup r-. 

veVi Iplivi 

The local version of m, applied to both primal and dual residuals, namely 


u-Ums\\v, < U||i?“||^. (A-Vl)-'/^ <C\\Rn\v, (A^+i)-'/" 


‘^msWy^i — ^ W-^H WV* V^Vi + 1 

motivates the local error indicator, defined as 

rif = mh^dmvrw::+i)-^ 


( 20 ) 

( 21 ) 

( 22 ) 

(23) 


where ||i?“||y.* and ||i?f ||y.* are defined in (fl^ and (1^ . respectively. Applying (1^ and (1^ to (ITOl) bounds 
the error in the goal function by 


g{u - Ums) < C^rjl. (24) 

i=l 

In summary, the goal-error over the global domain D is bounded by the product of energy errors of the primal 
and dual problems, which is in turn bounded by the sum of the indicators given by (1231) . modified by the 
partition of unity functions. This upper bounds suggests the adequacy of the indicators in reducing the error. 
The efficiency and a formal convergence analysis are however not addressed here. As in where a similar 
indicator is used for hp-refinement, this indicator displays similar behavior to the DWR-type indicator, as 
shown in the numerical experiments; however, it is more amenable to analysis. This indicator has the added 
advantage of reduced computational cost as compared to the DWR-type indicator described below, as both 
primal and dual problems are solved over the same discrete spaces, whereas for the DWR-type method, the 
dual problem has greater computational complexity than the primal. 
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DWR-type goal-oriented indicator 

The next error indicator is similar to DWR error indicator. For the primal problem solved in discrete space 
the DWR indicator is motived by the following residual equation. For 2 ; the solution to (O, it the 
solution to and Ums S Kff the solution to (|3]) 

g{u - Ums) = a{u - Ums, Z - Zos) = R'^{z - Zoff), (25) 

where Zoff in Voa is arbitrary and the global residual R'^{v) = K{x)Vums ■ Vi;. We let z*g 

be the component of ZoS spanned by the basis functions corresponding to the coarse neighborhood uji. 
Localizing (1^51) by the partition of unity functions Xi, 

N N N 

R-{z - Zoff) = ^ Rr = E (<nrich ” ^off) + E " Enrich) ■ (26) 

i—1 i—1 i—1 

As the exact solution z is unavailable, one generally instead replaces z by Zenrich, a discrete solution from 
a more enriched space than the primal, essentially neglecting the last term of (1^ . The function is 

the component of Zenrich spanned by the basis functions corresponding to the coarse neighborhood uii. In 
standard finite element methods the global residual is then solved elementwise and used as an indicator. By 
Galerkin orthogonality, then function Zoff may be taken as any function in Vqs but in practice is taken as 
the projection of the enriched dual solution into Voq. 

In this case, the dual problem is solved in the enriched space, called I4nrich- The space 14nrich is obtained 
by adding more basis functions to each coarse neighborhood. Recalling the construction of basis functions 
ipx,k in cni), the enriched space is constructed with more than li basis functions per coarse neighborhood, 
specifically 

'^i,k = for 1 < * < and 1 <k <li + ra, (27) 

where m basis functions are added for each Wi. The span of these basis is our I4nrich- Let Zenrich G Knrich 
be the solution for the dual problem in I4nrich, that is, Zenrich satisfies 

u(r^, Zenrich) — di.^), ^ L^nrich- (28) 


The DWR-type error estimator is defined as 


Vi = 


Ri ^Ti(Zenrich) (Zenrich)) ^ 


i = 1,2,--- ,iV. 


(29) 


In the above definition, Pi(zenrich) is the component of Zenrich spanned by the basis functions tpix, k = 
1,2, • • • ,li + m, corresponding to the coarse neighborhood uji. Moreover 7r(Pi(zenrich)) is the component of 
T’i(zenrich) Spanned by the basis functions '4’i,k, k = 1, 2, • • • ,li, in the offline space. Comparison with (l26l) 
yields an heuristic bound, modulo the error term created by replacing z by Zenrich- 

Each of the two indicators, given by (1231) . and respectively, (1291) . may be implemented in an adaptive 
framework to determine which coarse neighborhoods to enrich. The goal-oriented variant of the adaptive 
enrichment algorithm in Section [3.41 is now described. 


Goal-oriented adaptive enrichment algorithm 

Choose a fixed marking parameter 0 < 0 < 1. Choose also an initial offline space V^g by specifying a fixed 
number of basis functions for each coarse neighborhood, and this number is denoted by Ij. Then, generate a 
sequence of spaces and a sequence of multiscale solutions obtained by solving ([4]). Specifically, for 
each m = 1,2, • • •, perform the following calculations: 

Step 1: Find the multiscale solution in the current space. That is, find G such that 


a{uZ,v) = {f,v) foraWveV,^. 


(30) 
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Step 2: Find the multiscale dual solution in the current space or an enriched space. That is, find 
such that 

a(zZ,v) = (f,v) for all ^ e Kiu,i (31) 


where 


T/m _ 

*^dual 


off ’ 


enrich’ 


for JJ ^-based error estimator, 
for DWR-type error estimator. 


Step 3: Compute the local residual. For each coarse region uji, compute 




\\mv;A\Rnv:{xr:+ir\ 

|i?“(Pj(2;enrich) “ 7r(Pj (Zenrich))) |, 


for H ^-based error estimator, 
for DWR-type error estimator. 


R'iiv) = [ fv- [ ■ V^;, 

J UJi ^ OJi 

Ri{v)=f gv- I k(x)Vz^,-Vv, 

J U)i ^ OJi 

consistent with (HU and (1201) : and the norm is defined in (1131) and (1211) respectively. Next, re-enumerate 
the coarse neighborhoods so that the above local residuals rjf are arranged in decreasing order 77 ^ > 
77 I > • • • > 77 ^. That is, in the new enumeration, the coarse neighborhood wi has the largest residual 
77 ^. As in Remark l3.II the full sort of the estimators can be replaced by a partial sort for a marked set 
of quasi-optimal cardinality. 

Step 4: Find the coarse regions where enrichment is needed. Choose the smallest integer k such that 

N k 

(32) 

i=l i=l 

The coarse neighborhoods wi,a; 2 ,--- ,Wfc, are then enriched with additional basis functions. Alter¬ 
nately, a set of neighborhoods based on the binning strategy for a partial sort can be chosen as 
described in Step 3 of the Adaptive algorithm 13.41 

Step 5: Enrich the space. For each i = l,2,---,fc, add basis functions for the region uii according to the 
following rule. Let s be the smallest positive integer such that is large enough compared with 

Ajm+i. Then include the eigenfunctions • • • , in the construction of the basis functions. 

The resulting space is denoted as . Mathematically, the space is defined as 

= ^off + spanU^,^! 

where iliij = and 7 /;“’’°® = X]r=i j j = Z™ -|- 1 , • • • , -|- s, denote the new basis 

functions obtained by the eigenfunctions • • • , In addition, set = ^1" + ■s- 

In the next section, we demonstrate the efficiency of the goal-oriented adaptive algorithm defined above 
on a problem with high-contrast multiscale coefficients. The results are compared with the standard residual- 
based adaptive method defined in the previous section. We note both the increased efficiency in the reduction 
in goal-error, \g{u — Mrns)|, and the decreased reduction in the energy-norm error with the goal-oriented 
methods. These two observations suggest the method does what it was designed to do: focus the adaptive 
enrichment towards reduction in goal-error without resolving the solution where it has limited influence on 
the goal-error. We also note similarity in performance between the two indicators, both demonstrating errors 
with a similar observed rate of convergence. 
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5 Numerical Results 


In this section, we present two numerical examples for multiscale problems with high-contrast coefficients and 
compare the performance of the two indicators defined in the previous section. For our simulations, we take 
the domain = (0,1)^, and the inflow-outflow source term / = xki — where Ki = [0.1,0.2] x [0.8,0.9], 
and K 2 = [0.8,0.9] x [0.1,0.2]. We consider the problem of finding g{u) for u the solution to ((T|), namely 

—div(K(a;)VM) = /, in Z?, u = 0 on dD. (33) 


The goal functional 


9{u) 



u, 


is the average value of u on the outflow region K 2 . In practice, K 2 is the location of the wells, and it is 
important that the average pressure u on K 2 is accurate. The two examples differ by the high-contrast 
coefhcients k{x), shown in Figure O Shown on the left, ki features a high-conductivity channel crossing 
the domain separating the inflow and outflow; on the right, K 2 is a similar coefficient without the channel. 
These coefficents are visualized with the blue region indicating the value 1 and the red region indicating 
the contrasts. In each example, we consider two different contrast strengths, 10'’’ and 10®. We note the 
invariance in the relative performance of the indicators with respect to the contrast strengths. 



Figure 2: Left: the coefficient ki, corresponding to Figure El Right: the coefficient K 2 , corresponding to 
Figure m 

For the hrst example. Figure O shows a comparison in error reduction between the three indicators 
for coefficient ki. In the two figures on the right, we compare the logarithms of the energy norm errors 
g{u — Mins) against the number of unknowns for three types of adaptive enrichment algorithms; namely, 
the residual based method described in Section EH] (denoted in blue in Figure El), the H~^ residual-based 
goal-oriented method (denoted in red in Figure El, and the DWR type goal-oriented method (denoted in 
black in Figure E]) as described in Section [U From these results, we see the two types of goal-oriented 
methods behave similarly and outperform the standard adaptive method with an improved rate of goal-error 
reduction. We note a more stable decrease in error reduction for the goal-oriented residual-type method, 
but slightly improved, if less predictable error reduction for the DWR-type. This last observation is to be 
expected, as the DWR-type indicator does not account for the error created by using the enriched solution 
-^enrich in place of the exact dual solution 2 , as in (E51) . 

On the left of FigureEl we see the standard residual based method outperforms the goal-oriented methods 
for the energy norm error, ||m — MmsHy. This confirms that the goal-oriented methods are driving the adaptiv¬ 
ity toward a more efficient evaluation of the quantity of interest without expending additional computational 
effort resolving features of the solution with limited influence on the goal. 
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Figure 3: Lower contrast (10^). Top left: log ||u — UmsWv- Top right: log \g{u — Ums)\- 
Higher contrast (10®). Bottom left: log ||m — Ums\\v, Bottom right: log \g{u — Ums)\- 



DOF-DOF. 



DOF-DOF. 




Figure 4: Lower contrast (10^). Top left: log ||u — UmsWvi Top right: log \g{u — Ums)\- 
Higher contrast(10®), Bottom left: log ||u — UmsWvi Bottom right: log \g{u — Ums)\- 
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For the second example, we consider finding g{u) for u that satisfies (1551) . with k{x) given by K2 shown 
on the right of Figure[2l with no high-conductivity channel. As seen in FigurelH the results are qualitatively 
similar to the results of the first example. In summary, the plots on the right show the goal-error reduction for 
the higher and respectively lower contrast cases. The residual-type and DWR-type goal-oriented indicators 
achieve a better rate of error reduction than the standard adaptive method, with the DWR-type showing 
generally the lowest error, with the least-steady decrease. The plots on the right of Figure 0] show the 
reduction of energy error for the three indicators. As in the first example, the standard H~^ residual based 
adaptive method designed to reduce the energy error shows the best performance here, whereas the goal- 
oriented methods yield steady error reduction but are focused on localized error reduction in the region of 
the goal-functional, rather than across the entire domain. 

These results demonstrate the importance of goal-oriented adaptivity in GMsFEM, particularly in cases 
where the global domain is significantly larger than the region of infuence for the quantity of interest. In 
partiuclar, problems with many localized features only some of which significantly influence the quantity of 
interest will benefit from a goal-oriented adaptive strategy. 

6 Conclusion 

In this paper we defined two types of error indicators that can be used in an adaptive algorithm for multiscale 
problems with high-contrast coefficients. The goal-oriented adaptive algorithm fits within the framework of 
GMsFEM, and focuses the adaptivity on reducing the error in the quantity of interest, rather than in 
global norm. We first reviewed the general ideas of GMsEEM for high-contrast problems, then gave a 
detailed overview of the construction of multiscale basis functions, and a residual based adaptive algorithm 
designed to reduce the energy norm error. We stated the dual problem, then motivated and introduced 
two goal-oriented error indicators, and described their use in a a goal-oriented adaptive algorithm. Finally, 
we demonstrated the efficiency of the goal-oriented algorithm and estimators compared with the standard 
adaptive method introduced earlier. We found for both indicators, the goal-oriented method reduced the 
error in the goal-function at a better rate than the standard method. We also found the two indicators 
perform similarly, with some increase in error reduction seen in the DWR-type indicator, but at the cost of 
solving the dual problem in a more enriched space, increasing the computational complexity. The residual- 
based indicator on the other hand may be more amenable to convergence analysis, as may be investigated 
in future work. The current results indicate the goal-oriented strategy increases the efficiency of GMsFEM 
when a function of the solution rather than the solution in its entirety is of interest. 
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